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I. INTRODUCTION 

Spectacular creation of a Bose-Einstein condensate 
(BEC) in atomic traps [l| has opened unique possibil- 
ities for controlled studies of ultracold bosonic gases 0- 
One of the basic problems in physics of a BEC is how 
density perturbations propagate throughout a bosonic 
cloud. There has been a big interest in this problem 
and both an experiment 0] and theoretical studies 0, Q 
have been performed. These theoretical calculations have 
been done, however, for infinitesimally small perturba- 
tions, which dynamics hardly exhibits nonlinear effects 
leading to formation of shock waves. This paper is de- 
voted to theoretical studies of a shock wave formation in 
homogeneous and harmonically trapped BEC. 

Shock waves have been widely investigated in differ- 
ent physical systems GJ. For example, we find them in 
a gas bubble driven acoustically 0, in a photonic crys- 
tal and even in mathematical models of traffic flow 
@ . Of particular importance are studies of shock waves 
in classical compressible gases, where they appear when- 
ever large enough density perturbation propagates |^.ITo|. 
Contrary to the case of classical gases, shock waves in ul- 
tracold "quantum" atomic gases still require both theo- 
retical and experimental investigations even though some 
results have been already obtained. The propagation of 
steep density perturbations in a BEC has been exper- 
imentally investigated in the theory of shocks in 
Fermi (Tonks) gas has been proposed in [lj i an approx- 
imate theoretical descri ptio n of shock waves in a BEC 
has been worked out in |l3j |. It is also worth to mention 
that shocks have been studied in physical systems gov- 
erned by nonlinear Schrodinger equation |14lll5| . These 
studies can be related to investigations of shock waves 
in a BEC due to formal similarity between the nonlinear 
Schrodinger equation and mean-field equations of a BEC. 
Nevertheless, it is important to realize that dynamics of 
shock waves even in weakly interacting Bose gas, may 
differ from predictions of the mean-field approximation, 
as will be proposed in Sec. II. 

The mean-field description of a BEC is provided by 



the Gross-Pitaevskii equation [T^, which in the hydro- 
dynamical variables and dimensionless form |T^j reads 
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First equation is a continuity equation, while the sec- 
ond one is similar to the Euler equation from classi- 
cal hydrodynamics. Naturally, p is atomic density and 
v is a velocity field. The atomic density in JIJ satis- 
fies J dx p = 1, a constant g > is proportional to 
the strength of two-body repulsive interactions between 
atoms. The term cx d^^fp/ ^/p is called quantum pressure 
(QP). It is comparable to the gp term only if the spatial 
scale of density variations is of the order of the healing 
length £(p). Simple estimation on the basis of 0} shows 
that = l/y/Zj&. 

In the following wc assume that a BEC is initially in 
a state having a Gaussian-like density wave packet cre- 
ated by adiabatically focused laser beam on an atomic 
cloud. Alternatively, an atomic cloud can be cooled in a 
presence of a laser beam, as proposed in Q . We assume 
that the width of a laser beam is large compared to the 
system's healing length, so that the density perturbation 
is broad and the QP term, at least initially, is negligible. 
We are interested in dynamics of such a laser-induced 
perturbation after abrupt laser turn off. We show that 
it splits into two wave packets propagating in opposite 
directions. These wave packets undergo self-steepening 
dynamics. Finally, two shock wave fronts are formed. 



II. ID HOMOGENEOUS SYSTEMS 

First we consider the case of a one-dimensional (ID) 
BEC in a periodic box having boundaries at a; = ±i. 
There is initially a laser beam producing the potential 
V = uq exp(— x 2 /2a 2 ) with a <C I and uq being controlled 
by intensity and detuning of a laser beam. The initial 
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density of atoms takes form 

p(a;,0) = po + 2w^ /2 ' 2 



1.5x10 



(2) 



where po = 1/(21 + 2r\\/2/Ea) provides normalization of 
iJSJ. We assume that the system is in the ground state 
which implies v(x, 0) = 0, and that the QP term is neg- 
ligible: a ^> £,(po). It means that we are in the so-called 
Thomas-Fermi (TF) regime Under such conditions 
2f?Po = -u Q /g. 

Let us first inspect what happens in a low amplitude 
limit. The analysis can be done in the framework of a Bo- 
goliubov approach or within a hydrodynamical approach. 
The first possibility was explored for ultracold fermions 
in ^3 ■ Here we follow the second one to give the reader 
more insight into the problem. 

We express p(x,t) as po + Sp(x,t) with \Sp(x,t)/po\ -C 
1. Having in mind that the velocity field v(x, t) is of the 
same order as Sp(x,t), a linearization of Q without the 
QP term gives 
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where cq = ^gpo is a speed of sound wave propagat- 
ing on density po [l8j • Straightforward calculation deter- 
mines the solution of © that satisfies initial conditions 
5p(x, 0) = 2rjp exp(— x 2 /2a 2 ) and v(x, 0) = 0: 

cj)(x,t) = [ dke- k2a2/2 g[sin(kx + ) - sm(kx-)]/kcf), 

v27T J 

with x± = x =p cot. Using the explicit form of 4 > {x 1 1) and 
© we find 

Sp(x, t) = VPo [ e -(- c °*) 2 / 2ff2 + e -(x+c»t) 2 /2^ . (4) 

Therefore, the initial perturbation splits into two pieces 
traveling in opposite directions with a speed of sound 
during free time evolution (Fig. Qi-b). Such a process 
is impossible in a system which dynamics is governed 
by a single-particle Schrodinger equation. It is a nice 
interaction- induced phenomenon in a BEC. Interestingly 
a similar effect can be observed in Fermi (Tonks) clouds 
[l2|. It is also easy to see from (|2J| that the splitting 
phenomenon happens in exactly the same way regardless 
of the form of an initial shape. 

Let us also remark that, since we consider particles in a 
periodic box, all the equations should be strictly speaking 
written in an appropriate Fourier basis. For a <C / and 
c$t < I the discrepancy between obtained results and 
exact ones is hardly visible. 

Even though the solution (QJ was derived for small per- 
turbations (|?7| <C 1) it can be found numerically that the 
splitting of an initial wave packet into two oppositely 
moving pieces takes place for |?7| ~ 1 as well. Quantita- 
tive comparison of Q to exact numerical results reveals 




FIG. 1: Density of atoms during free time evolution. Dotted 
line - an initial density profile, dashed line - the prediction 
@, solid line - a numerical solution of Q with the QP term. 
Plot (a) [(b)]: rj = 0.0175 and t = 21 [t = 42]. Plot (c) [(d)]: 
the right-moving wave packet for r\ — 0.0175 \q = —0.0175] 
at t — 63. Other parameters: g = 7.5 • 10 3 , a — 6.5, I = 350. 



discrepancies growing in the course of time evolution for 
any \rj\. An exact solution for rj > (?/ < 0) moves faster 
(slower) and changes its shape so that the front (back) 
impulse's edge self-steepens - Fig. ^>d. 

To explain that behavior we have found, following 
the exact traveling wave solution of Eqs. without the 
QP term: 



P = f(x~ (a Q ±3y/gp) t) 



ao ± 2y/gp, (5) 



where / (an arbitrary function) and ao (a constant) can 
be determined from initial conditions. The sign ± deter- 
mines a direction of propagation. The result (|SJ is valid 
as long as a spatial scale of density perturbations is larger 
than the healing length £. Solutions of the form JSJ were 
previously worked out in various physical systems, see e. 

g. dEiiia. 

Nonlinearity of Q prevents usage of the superposition 
principle to built a solution, that contains both left and 
right-moving perturbations, out of JSJ. Fortunately, the 
problem can be simplified, since after separation pertur- 
bations move independently. 
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We propose to consider only the right-moving 
(r) part and to take as an initial density profile 

p r (x, 0) = po + poV e ~ x ~^ 2cr2 (compare ©)■ Such an ap- 
proximation is exact for small perturbations and turns 
out to be good for large ones - see discussion of numerical 
simulations presented below. The constant ao is found to 
be —2y/gpo from a requirement that initially the velocity 
field far from the perturbation is zero: v(x — > ±1, 0) = 0. 
This assumption leads to the solution 

p r (x,t) = p +p r]exp(-[x-y^(-2 v ^+3 v ^)t\ 2 /2a 2 ). 

(6) 

Although ® is in an implicit form, all the quantities 
characterizing a shock wave formation can be analytically 
extracted from it. 

Both amplitude and velocity of impulse's maximum are 
constant during time evolution. The amplitude equals 
p (l + if), while the velocity becomes V(i]) = ^/gpoA(r/) 
where 



A(rj) = 3^1 + ^-2. 



(7) 



For small perturbations we get V(rf) — ^/<?po(l + 'irj/2). 
Taking r\ = we find a well known expression for a sound 
velocity: ^/gpo [),&]■ These results show explicitly that 
bright perturbations (77 > 0) move faster than dark ones 
(r) < 0). A typical measurement of a sound velocity relies 
on observation of propagation of density perturbations on 
a cloud . If such perturbations are not small enough 
our results indicate that the r/ dependent term in V{rf) 
has to be considered in the experimental determination 
of a sound velocity. 

For clarity of discussion we will concentrate now on 
bright perturbations. A speed of impulse's maximum is 
bigger than a speed of its tails. As a result an im- 
pulse self-steepens in the direction of propagation so that 
formation of a shock wave front takes place. Finally, 
d x p r (x,t s ) = — oo at x = x s leading to the well known 
wave breaking phenomenon [6j. 

Time required for such a process (t s ) can be estimated 
as follows: the impulse's half-width 2a) is a difference 
in distance traveled by lower and upper impulse's parts 
until the shock appears: (V(?7) — V(0))t s w 2a. It gives 
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Quantitatively time of shock creation and shocks' posi- 
tion can be determined from a set of equations [j| 
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The first one says that d x p r is infinite, while the second 
one assures uniqueness of the p r {x,t s ) function. Assum- 
ing r] > we find that 



x iPr) = \/5(-2\/Po+ 3 \/ft : ) t+y/-2a 2 \n[(p r /p Q - l)/rj\], 

(10) 
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FIG. 2: Schematic plot of functions: —In 7 (dashed line) and 
(l+»77)/ (2 +777) (solid line) vs. 7. Dotted lines define bounds 
leading to inequalities 1131 . 



which substituted into © leads to 

y/2a y/p s + po 



ts = 



3V5 Ps - Po 

where p s , the density at which d x p r — —00, satisfies 
- In [(Ps/po - 1)/??] = ■ 



(11) 
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It seems to be impossible to solve analytically that equa- 
tion. Fortunately, its numerical treatment is straightfor- 
ward. It is possible, however, to prove analytically that 
the solution of l|12[) exists and satisfies the relation 



po + pove B{11) <p s <Po + Porie- 1 ' 2 , 



(13) 



where 
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For 77 —>■ lower and upper bounds coincide. One can 
derive inequalities l|13f) by means of a graphical method 

- Fig. |21 Let us introduce a variable < 7 < 1 such 
that p s = po(l + 777). Eq. I|12l) now takes the form: 

— In 7 = (1 + 777) f(2 + 777). The RHS of this equation in- 
creases monotonically with 7 from 1/2 to (1 + 77) / (2 + 77). 
Conversely, the LHS decreases monotonically with 7 from 
+00 to 0. Therefore there must be 7 at which the RHS 
equals the LHS. For 7 > exp(— 1/2) the LHS becomes 
smaller than 1/2. It means that 7 < exp(— 1/2) in the 
considered problem. Therefore, the RHS takes the high- 
est value for 7 = cxp(— 1/2), so it is bounded from above 
by B(rj). As a result, we arrive at a conclusion that the 
solution exists and 1/2 < —In 7 < B(rf) leading to in- 
equalities i|13|) . 

On the basis of l|ll|) and i|13[l we have found bounds 
for time of shock formation 



ts < ^ ^V2- 
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In the 77 — > limit both bounds equal 2 ere 1 / 2 /3rjy/gpo. 
From inequalities (|14|1 we see that the shock appears if 
a wave packet evolves long enough. This observation to- 
gether with the explicit expression for a speed of propaga- 
tion of density maximum is a central result of this paper. 
Considerations ©-(01 fail to predict such behavior since 
they are valid in the 77 — ^ limit where t s — > 00. As 
mentioned above these calculations have been done for 
bright wave packets (77 > 0). Extension to dark impulses 
is straightforward. We would like also to point out that 
the solution JSJl is uniquely defined in a whole space up 
to the time of shock creation, i.e. for t < t s - see Fig. 0|d 
and Ref. . Later on, the solution © in a shock wave 
front region becomes multi- valued, compare 0, El E3 • 
Approximating the position of a shock wave front by 
position of impulse's maximum at t = t s : x s = V(rj)t s , 
we get 



x s > 



377 

377 



v / 2 + 7;e- 1 /M(?7). 



(15) 



As rj — > we have x s — > 2ere 1 / 2 /3rj. 

Let us compare numerical simulations of hydrodynam- 
ical equations with the QP term to above derived an- 
alytical predictions. We must be aware that analytical 
calculations have been performed without taking the QP 
term into account. Initially a spatial scale of density 
variations is roughly a, which is much greater than the 
healing length £. As a result the QP term is negligi- 
ble at the beginning of time evolution. As shocks form 
up density changes occur on smaller and smaller length 
scales being finally of the order of £. Therefore, before 
appearance of shocks quantum pressure has to start to 
play a role leading to discrepancy between analytical pre- 
dictions and numerical simulations. 

For numerical purposes, we have fixed the dimension- 
less coupling constant g to 7.5 • 10 3 . In a typical 3D con- 
figuration such a value corresponds to roughly 10 5 ~ 10 6 
atoms of either 87 Rb or 23 Na. In the present calcula- 
tions such a value of g guarantees that we are well in the 
Thomas-Fermi regime, which is important for discussed 
phenomena. The same qualitative results were obtained 
for g in the range 2 • 10 3 — 5 • 10 4 . 

We have performed numerically time evolutions and 
have observed density wave packets until they reached 
box boundaries being at x = ±1. These results can be 
divided into two groups: the case of x s > I and that of 
x s < I. In the first situation we do not observe shocks, 
although a shock-like deformation of an impulse's shape 
can be clearly visible. The agreement between analytical 
solution Jfjjl and full numerical simulation is very good - 
Fig. [3J When x s < I shocks appear. A typical situation 
is depicted in Fig. H^-c. 

First of all, as a relative height of a wave packet grows, 
amplitude of separated density perturbations becomes a 
little bit smaller than half-amplitude of an initial pertur- 
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FIG. 3: Density of atoms during free time evolution. Dotted 
line - an initial density of atoms, dashed line - the analytical 
solution @ , solid line - a numerical solution of Eq. Q with 
the QP term. Plot (a) [(b)] is for t = 24 [t = 51]. Other 
parameters: 77 = 0.06, t„ = 60.4, x s = 253.8, g = 7.5 • 10 3 , 
a = 12.5, I = 250. 



bation. It leads to a small discrepancy between analytical 
predictions and numerics (Fig. UJi). 

Second, a slope of the impulse does not become infi- 
nite at t = t s , but significantly increases. To express it 
quantitatively we have defined a relative slope by 



S(t) 



max(- 



dp(x,t) * 
dx > 



-1/2 



(16) 



where the denominator is max(— for a single wave 
packet which density is p + 7/p exp(— x 2 /2cr 2 ). After 
the splitting 3»1 and than grows. As depicted in Fig. 
0J; the relative slope can be as large as few tens showing 
that a shock formation occurs on a time scale t ~ t s . 

Third, for t > t s density oscillations in the mean-field 
approach with the QP term show up instead of a true 
shock wave front. As shown in Fig. for t > t s the 
analytical solution JHJ still fits to numerics, but in a back 
edge of an impulse only. In the front part the solution 
© becomes multi-valued, as expected from the theory 
of nonlinear equations 0, . We have checked that the 
appearance of density oscillations is a result of coming 
into play of the QP term. It is not an artifact of inability 
of our numerical methods to evolve a truly shock wave 
front if such would appear [2(]]. Qualitatively the same 
results were obtained by other means in [13|, [Tj, [l5| . 

Similar density oscillations have been recently observed 
in a paper that describes shock wave formation in a Fermi 
(Tonks) cloud 01 ■ We have shown there that they are 
an artifact of the usage of the mean-field approximation. 
Indeed, they were absent in exact many-body calcula- 
tions. In a Fermi (Tonks) cloud there were three stages 
of shock dynamics: a shock wave formation, propagation 
of a shock-like impulse roughly without change of shape, 
and impulse's explosion leading to broadening of density 
profile. The mean-field approach has reproduced only 
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FIG. 4: Plots (a) and (b): density profile of the right-moving 
wave packet. Dashed line - the analytical solution ©, solid 
line - a solution of Q with the QP term, thick solid line in 
(b) shows multi- valued part of the solution JSJ . Plot (a) [(b)] 
is for t = 18.9 [t = 27]. Plot (c): relative slope GJJ of tne 
right-moving wave packet. Dotted line indicates time of shock 
formation l|ll|l - H12|l . Parameters: rj = 0.2, a = 12.5, t s = 19, 
x s = 93.5, g = 7.5 ■ 10 3 , I = 250, £(2 ■ 10" 3 ) = 0.18. 



the first stage. The second stage was absent - density 
oscillations have appeared in the same way as here. As a 
result, we speculate, that it is rather doubtful that pre- 
dictions of the Gross-Pitaevskii equation at the front edge 
are reliable after appearance of density modulations. 

It is difficult to make any definite statements about 
applicability of the Gross-Pitaevskii equation for descrip- 
tion of shocks without the insight into a many-body the- 
ory of interacting bosons. To understand the possible 
source of problems we recall that the N-body bosonic 
wave function has a general form: 

\/\q4>(xi) ■ ■ ■ 4>(xn) + an orthogonal part. 

The condensate fraction Ao € (0, 1] is divided by N the 
highest eigenvalue of a single particle density matrix. 
4>{x) is the eigenvector of that matrix to eigenvalue Ao- 
In the mean-field approach we assume Ao ~ 1 . Neverthe- 
less, Ao can change during time evolution. In principle, it 
can be noticeably lower than one after a shock creation. 
It means that contribution of states orthogonal to a con- 
densate mode, neglected in a mean- field treatment, might 
be significant. These states may smooth density changes 
by filling the oscillatory region with atoms depleted from 
a condensate. Similar situation one encounters for dark 
soliton in a BEC [2^, namely a soliton's notch, which 
size is also of the order of the healing length, fills with 
depleted atoms. 



III. ID HARMONICALLY TRAPPED SYSTEMS 

Now we would like to find out what happens if an ex- 
ternal harmonic trapping potential is imposed on a ID 
bosonic cloud. Atoms are initially placed in both har- 
monic and laser potentials 01 

x 2 

V = — + u exp(-x 2 /2<j 2 ). 
The atomic density, in the TF limit, takes form 

p( x ,0) = H (|.7) 

g 
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^0_ e -x 2 /2o 
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for |x| < Rtf — \/2p an d zer0 f° r M > Rtf- Similarly 
as before, relative amplitude 2rj of an initial perturbation 
is defined by relation 2rjp/g = —uo/g. We have restricted 
our numerical simulations to bright perturbations in this 
section: r\ > 0. We have also assumed that the width of 
a laser beam is much smaller than spatial extent of the 
cloud. Under these conditions the chemical potential p, 
found from the normalization of l|17(l . equals 

p(u ) = [(1 + u V2^a/g)3g/4V2] 2 / 3 . (18) 

Once again, we are interested in dynamics of Gaussian- 
like wave packets created after an abrupt laser turn off. 

First of all, a splitting of an initial perturbation into 
two propagating outward pieces takes place regardless of 
an initial position of a density perturbation. Secondly, 
there appears self-steepening of the bright (dark) im- 
pulse's front (back) edge - see Fig. [SJi. Finally, den- 
sity oscillations appear showing that an impulse enters a 
shock regime. 

We were unable to find exact solutions of Eqs. (JTJ in 
a trapped case even when the QP term was absent. In- 
stead of it, we can provide a simple semi-analytical pre- 
diction for a speed of density perturbations on the cloud. 
To this aim we add x 2 /2g term to the density of atoms 
p{x,t). Such a transformation removes a harmonic trap- 
ping term in (|17|l and propagation of impulses is qualita- 
tively similar to the homogeneous case carefully discussed 
above (Fig. \Sjp). Depending on amplitude of outward 
moving perturbations, in the p(x, t) + x 2 /2g picture, im- 
pulses' amplitude can change a little during propagation. 
Therefore, it is not easy to find out phenomcnologically 
a simple rule that governs motion of Gaussian-like wave- 
packets in a harmonically trapped system. 

We have tried different semi-analytical formulas for the 
prediction of a position of a wave packet's maximum, and 
found that the following approach is the most accurate. 
In a homogeneous case a speed of a density perturbation 
was greater than a local sound velocity by the factor A{rj) 
0, with r] being a relative impulse's amplitude ©■ We 
have assumed that a similar law holds also now. It leads 
to the following approximate equation for the position of 
a maximum of p{x, t) + x 2 /2g 



dx r 



dt 



A(rj)^/p-x 2 n /2. 



(19) 
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FIG. 5: Plot (a): density of atoms at different time moments. 
Plot (b): the same as in the plot (a) but with the x 2 /2g term 
added to the density profile in order to remove a trivial posi- 
tion dependence of atomic density due to external harmonic 
trapping potential - see Eq. I17H . Dotted line - 1 = 0, dashed 
line - t — 0.3, dashed-dotted - t = 0.6, solid line - t = 0.8. 
Both t and x are expressed in harmonic oscillator units [l7|. 
Other parameters: u = -90, rj = 0.184, g = 7.5-10 3 , a = 1.4, 
Rtf ~ 22. 



Integrating l|19|) we have found 



x m (t,r)) = -Rtf sin 



V Rtf 



(20) 



In the limit of r\ — ► the zero-amplitude result of p| 
is recovered. We have checked prediction (|20[l for 77 as 
large as 0.4 and g — 7.5 • 10 3 . Notice that it corresponds 
to relative amplitude of an initial perturbation equal to 
80%! Numerical results are depicted in Fig. E^-c. As 
easily noticed, the agreement between numerics and the 
formula l|20[) is satisfactory. We have observed that as 
x m — > Rtf, Eq. I20f) underestimates slightly a position 
of impulse's maximum for small 77 and overestimates ex- 
act result for large 77 j^. For 77 w 0.18 the agreement 
between numerics and semi-analytics is excellent - Fig. 
EJj. Despite differences far away from a trap center it is 
important to notice that the formula l|20|l correctly works 
near a trap center, where density is quite homogeneous. 
Indeed, Eq. I|20(l reproduces nicely numerics at least for 
x m {t,rj) < Rtf/2, V € (0.02,0.4), and other parameters 
as in the caption of Fig. |SJ 

On the basis of l|20|l we conclude that near a trap cen- 
ter, where Rt f sm(A(j]) ^/Jlt / Rt f) ~ A(i])y/Jlt, the P°" 
sition of impulse's maximum changes linearly in time and 
its speed of propagation equals A{rj)yfjl. This observa- 
tion is in qualitative agreement with the experiment per- 
formed in a quasi-one-dimensional configuration |3( • The 
most important message from the formula (|20|l is that the 
speed of finite size wave packets can differ significantly 
from the zero-amplitude result by the factor A{rf). For 
a realistic value of rj — 0.3 a speed of a wave packet's 
maximum is greater from the zero-amplitude prediction 




FIG. 6: Plots (a)-(c): a position of impulse's maximum vs. 
time. Solid line - the semi-analytical prediction 1201 , dashed 
line - numerics on the basis of Q with the QP term, dot- 
ted line - time of shock formation extracted from numerics 
(see text). Dashed-dotted line on the plot (b): x m (t,0). Pa- 
rameters (1*0,77) for plots (a)-(c): (-30,0.06), (-90,0.184), 
( — 180,0.39), respectively. Plot (d): time of shock formation 
vs. relative amplitude of an impulse. Solid line comes from 
numerics, while dotted line is an approximation l|2H . Both t 



and x are expressed in harmonic oscillator units 



Other 



parameters: g — 7.5 • 10 , a = 1.4, Rtf ~ 22. 



by more than 40%! We would like to stress, however, 
that these estimations are obtained for idealized case 
of a one-dimensional BEC, while the real experiments 
are typically performed in quasi-one-dimensional setups. 
Further studies are needless for clarification of whether 
similarly strong dependence of the speed of finite size 
wave packets on their amplitude survives in a quasi-one- 
dimensional configurations similar to that of |3j. The 
propagation of tiny density perturbations in harmoni- 
cally trapped quasi-one-dimensional setups was discussed 
in 0i 13 • These papers show that the tight confinement 
in cigar-shaped configurations changes the sound velocity 
expected for one-dimensional harmonically trapped sys- 
tems by a factor 1 / \[2. If the same property would hold 
for finite density perturbations, the application of our re- 
sults from this section to quasi-one-dimensional systems 
would be straightforward. 

From our simulations we have extracted time of shock 
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formation, defined now by a time moment when density 
oscillations in a shock front have amplitude large enough 
to lower density p(x) + x 2 /2g below the background value 
fi/g. We have found that shock-like behavior shows up 
for any r\ € (0.06, 0.45) and other parameters as in the 
caption of Fig. H3 The shock behavior of a large wave 
packet takes place in a relatively homogeneous part of 
the cloud where a sound velocity is ~ y/Jl. Replacing in 
the formula JSj ^gpo by y/JI the following estimation for 
time of shock creation can be derived 



s 3(yr+^-i)VM' 

As depicted in Fig. [BJi, this semi-analytical formula 
nicely reproduces time of a shock formation for large 
enough 77, and parameters of Fig. despite its approx- 
imate character. The same good agreement is found for 
g = 10 4 and g = 1.25 • 10 4 . For g lower than 7.5 • 10 3 , 
the case presented on the plot, discrepancy between es- 
timation and exact calculation grows up, e.g. i(2T|) 
overestimates time of shock formation at r\ = 0.4 and 
g = 2.5 ■ 10 3 by about 10%. 

IV. CONCLUSIONS 

We have considered propagation of broad density per- 
turbations in a BEC. Our theoretical approach fully ac- 
counts for nonlinear effects in sound propagation so it is 
well beyond the common zero-amplitude approximation 
0, , and allows for description of shock formation. As 
a result, quantitative predictions concerning propagation 
of large density perturbations in a homogeneous system 
were given. Extension of these results to harmonically 
trapped problems resulted in a semi-analytical descrip- 
tion being in a satisfactory agreement with numerics. 

We have focused on ID systems recently realized ex- 
perimentally |24| , but we expect that our work may pro- 
vide a theoretical basis for investigations of quasi-one- 
dimensional setups widely studied in laboratories. The 
easiest for experimental inspection are predictions for 
harmonically trapped condensates. We expect that also 
findings for homogeneous bosonic clouds can be verified 
with the help of, e. g. atomic wave guides pEj . 

It is worth to point out that shock-like structures have 
been already studied experimentally in a two component 
BEC H3> where time evolution of dark narrow density 
perturbations was observed. The initial size of a wave 
packet in jll| was of the order of the healing length, 
which made observation of impulse's self-steepening, a 
clear signature that shock formation takes place, very dif- 
ficult. We expect that the experimental setup proposed 
by us should give more convincing results. In fact, it 
seems that a slight change in a setup of together with 
a little bit higher measurement accuracy should verify 
our predictions at least qualitatively. 

From the theoretical side we would like to note that 
formation of shock-like solutions in a harmonic trap has 



been noticed in |4|. Our considerations provide an an- 
alytical insight into their numerics. It is also worth to 
mention, that it is possible to find an approximate an- 
alytical solution of hydrodynamical equations with 
the quantum pressure term. That solution also exhibits 
shock behavior, and its properties have been discussed in 

PHI HI 

It is also interesting to discuss differences between the 
traveling wave solutions presented in this paper and re- 
cently investigated, both theoretically and experimen- 
tally, solitonic solutions in a BEC |2fj. Our paper de- 
scribes density perturbations which change their shape 
during time evolution, while soliton-like de nsity wave 
packets propagate in dispersionless manner hj. The 
main difference between time evolution of wave packets 
considered here and solitonic ones results mainly from 
the difference in their width. The width of a soliton is of 
the order of the healing length, which makes the quan- 
tum pressure term in Eq. Q considerable. This allows 
for finding a solution that propagate without changes of 
shape. Our solution, Eq. JSJ, describes a perturbation 
that is much broader than the healing length. In this case 
the quantum pressure term is negligible and instead of 
solitonic solutions shock waves show up. It is also worth 
to realize that the initial shape of a shock-like solution, 
determined by the function f(x) in Eq. (JSJ), can be ar- 
bitrary as long as f(x) changes significantly on length 
scales much larger than the healing l eng th. Conversely, 
solitons have precisely defined shape [H| . 

One of the main open questions is how the shock-like 
impulse's front evolves in time. We have pointed out that 
the mean-field approach gives rather doubtful answer and 
suggested a possible explanation. Further investigations 
to resolve that issue will be a subject of future stud- 
ies. Unfortunately, contrary to the case of Fermi (Tonks) 
gases |12|. it is very difficult to solve a many-body prob- 
lem exactly. 

Finally, it is worth to point out that both a splitting 
of an initial perturbation and a shock wave formation 
take place in classical hydrodynamics 0]. It suggests 
that qualitatively the same phenomena can be observed 
in a cold bosonic (fermionic) gas above the condensation 
(Fermi) temperature. The crossover from a classical non- 
degenerate to a quantum degenerate regime can be an in- 
teresting subject for future experimental and theoretical 
investigations. In particular, differences in properties of 
cold bosons and fermions can be very exciting. Indeed, 
e. g. a sound velocity of ID ultracold bosons (fermions) 
is ~ y/p (~ p), while in the non-degenerate regime it has 
to have the same dependence on density for bosons and 
fermions. 
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